Synchronization of networks with variable local properties 
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We study the synchronization transition of Kuramoto oscillators in scale-free networks that are characterized 
by tunable local properties. Specifically, we perform a detailed finite size scaling analysis and inspect how the 
critical properties of the dynamics change when the clustering coefficient and the average shortest path length 
are varied. The results show that the onset of synchronization does depend on these properties, though the 
dependence is smooth. On the contrary, the appearance of complete synchronization is radically affected by the 
structure of the networks. Our study highlights the need of exploring the whole phase diagram and not only the 
stability of the fully synchronized state, where most studies have been done up to now. 

PACS numbers: 05.45.Xt, 89.75.Fb 



INTRODUCTION 

Emergent collective phenomena have been studied since 
long time ago. These phenomena arise in many fields of sci- 
ence, ranging from natural to social and artificial systems. 
They are characterized, among other features, by the collec- 
tive behavior of many interacting units that show a pattern 
hard to predict from the individual behavior of the system 
constituents. Several seminal models of statistical physics 
and non-linear dynamics have been scrutinized as paradigms 
of self-organization, emergence and cooperation between the 
units forming the system. In particular, synchronization phe- 
nomena constitute one of the most striking examples because 
of the many systems showing synchronization patterns in their 
behavior (Hill. 

One of the most celebrated synchronization models is due 
to Kuramoto 0, 0, who analyzed a model of phase oscil- 
lators coupled through a function (sine) of their phase dif- 
ferences. This model owes most of its success to the plenty 
of analytical insights that one can get through the mean-field 
approximation originally proposed by Kuramoto. In this ap- 
proach (KM), the nodes of an all to all, i.e. globally, coupled 
network, are considered to be oscillators with an intrinsic fre- 
quency and their phases evolve in time in such a way that if 
the coupling between them is larger than a critical threshold, 
the whole system gets locked in phase and attains complete 
synchronization. 

However, it has been recently discovered that real sys- 
tems do not show a homogeneous pattern of interconnections 
among their parts. That is, the underlying structure is not com- 
patible with the original assumption of the KM. It is not even 
well described by random patterns of interconnections in the 
vast majority of systems. Therefore, the mean-field approach 
requires of several constraints that are not usually fulfilled in 
real systems. Natural, social and technological systems show 
complex patterns of connectivity that characterize seemingly 
diverse social 01, biological 0,(3] and technological systems 
They exhibit common features that can be captured 
using the tools of graph theory or in more recent terms, net- 
work modeling ElElEl- 



It turns out that many real networks are well described 
by the so-called scale-free (SF) networks. Their main fea- 
ture is that the probability that a given node has k connec- 
tions to other nodes follows a power-law Pk ~ k^ 1 , with 
2 < 7 < 3 in most cases fill [Till . The study of processes 
taking place on top of these networks has led to reconsider 
classical results obtained for regular lattices or random graphs 
due to the radical changes of the system's dynamics when 
the heterogeneity of complex networks can not be neglected 

oehqgieiiii. 

It is then natural to investigate how synchronization phe- 
nomena in real systems are affected by the complex topologi- 
cal patterns of interaction. This is not an easy task, as one has 
to deal with two sources of complexity, the nonlinear char- 
acter of the dynamics and highly non trivial complex struc- 
tures. In recent years, scientists have addressed the problem of 
synchronization capitalizing on the Master Stability Function 
(MSF) formalism 11811 which allows to study the stability of 
the fully synchronized state (T^ E3 E3] - While the 

MSF approach is useful to get a first insight into what is going 
on in the system as far as the stability of the synchronized state 
is concerned, it tells nothing about how synchronization is at- 
tained and whether or not the system under study exhibits a 
critical point similar to the original KM. To this end, one must 
rely on numerical calculations and explore the entire phase 
diagram. Surprisingly, there are only a few works that have 
dealt with the study of the who le sy nchronization dynamics in 
specific scenarios 1251 12(1 1271 l28l |2jJ 13(11 as compared with 
those where the MSF is used, given that the onset of synchro- 
nization is readier in its behavioral repertoire than the state of 
complete synchronization. 

In this paper, we take a further step in the detailed charac- 
terization of the phase diagram and specifically, in the descrip- 
tion of the dynamical behavior at the onset of synchronization 
in SF networks. By performing a standard finite size scaling 
analysis, we show that the local topology affects the critical 
properties of the dynamics, though it is less pronounced than 
what one may expect a priori. We capitalize on a network 
model that keeps the power-law exponent fixed while varying 
the clustering coefficient and the average path length. In what 



2 



follows, we describe the topological and dynamical model and 
discuss the results from a global perspective. Finally, in the 
last section, we state our conclusions. 



NETWORK MODEL AND DYNAMICS 

We implement a network model in which the graph is 
grown at each time step by linking preferentially new nodes to 
already existing nodes in the same way as in the Barabasi and 
Albert (BA) model l3lll . The only difference is that the nodes 
are assumed to have a fitness that characterizes their affinities 
. In this way, by tuning a single parameter \i, one can go 
from the BA limit down to a network in which several net- 
work properties vary as a function of /i. On the other hand, as 
the linking mechanism is still the BA preferential attachment 
rule, the exponent 7 of the power-law degree distribution is 
the same (i.e., 7 = 3) regardless of the value of [i. Roughly 
speaking, the model mimics the situation in which new nodes 
are attached to an existing core or network but without having 
knowledge of the whole topology. 

The recipe is then as follows J32I1 . i) Initially, there is a 
small, fully connected, core of too nodes. Assign to each of 
these mo nodes a random affinity taken from a probability 
distribution. In this work, we have used for simplicity a uni- 
form distribution between (0, 1). ii) At each time step, a new 
node j with a random affinity dj is introduced and to links are 
established with nodes already present in the network follow- 
ing the rule 

n(fcO - h , (i) 

where the set Y contains all nodes that verify the condition 
di — fj, < dj < a,i + fi, being fx g (0, 1) a parameter that 
controls the affinity tolerance of the nodes. Finally, Hi) repeat 
step (ii) t times such that the final size of the network be N = 
mo + 1. 

In the above model, when \i is close enough to 1, the BA 
model is recovered. When it is decreased from 1, the values of 
some magnitudes such as the clustering coefficient ((c)) and 
the average path length ((£}) grows with respect to the BA 
limit I32I1 . In Fig.^ we have represented how these properties 
vary as a function of the parameter [i. Note that the larger vari- 
ations correspond to the clustering coefficient (a factor greater 
than 4 as compared to a factor close to 2 for (L)) and that it 
is the first property that deviates from the BA limit. This ten- 
dency holds up to very small values of fx, where (L) raises at a 
higher rate than (c) (not shown in Fig.Q. More important for 
our purposes is the region of 0.4 < fj, < 1. For these values 
of \x one observes that (L) remains constant while (c) starts to 
grow as soon as it moves away from the BA limit(/i = 1). This 
allows to decouple the effects of both magnitudes on what we 
are going to study. As we shall latter see, the structural clus- 
tering plays a major role in the synchronization of Kuramoto 
oscillators, as does in other dynamical processes 13311 . 
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FIG. 1 : (color online) Top: Values of the clustering coefficient rela- 
tives to those of the BA model against the parameter u. Bottom: The 
average cluster length for the network generated relatives to the BA 
value as a function of [i. All the networks are made up of N = 1000 
and have an average degree (k) = 6. 



The dynamic ingredient of the model is given by the collec- 
tive behavior that arises when the nodes are considered to be 
phase oscillators that follow the Kuramoto model. In this for- 
malism, the population of N interconnected units are coupled 
phase oscillators where the phase of the i-th unit, denoted by 
6i(t), evolves in time according to 

d0- 

=Wi + J2 A ij A iJ sin ( e 3 -h) * = 1,...,JV (2) 

3 

where u)\ stands for its natural frequency, Ay = A |34| is the 
coupling strength between units and Ay is the connectivity 
matrix (j4y = 1 if i is linked to j and otherwise). Note 
that in the original Kuramoto model mean-field interactions 
were assumed which leads to Ay = K/N\/i,j, for the all-to- 
all architecture. On the other hand, the model can be solved 
in terms of an order parameter r that measures the extent of 
synchronization in a system of N oscillators as: 

re» = l£y* (3) 
j=i 

where ^ represents an average phase of the system. The pa- 
rameter r takes values < r < 1, being r = the value of 
the incoherent solution and r = 1 the value for total synchro- 
nization. 
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FIG. 2: (color online) Order parameter rasa function of A for dif- 
ferent values of \a as indicated. The network parameters are those of 
Fig.0 

RESULTS 

In order to inspect how the dynamics of the N oscillators 
depends on the underlying topology, we have performed ex- 
tensive numerical simulations of the model. Starting from 
A = 0, we increase at small intervals its value. The natural 
frequencies and the initial values of 0i are randomly drawn 
from a uniform distribution in the interval (—1/2, 1/2) and 
(— 7r, 7r), respectively. Then, we integrate the equations of 
motion Eq. using a 4 th order Runge-Kutta method over 
a sufficiently large period of time to ensure that the system 
reaches the stationary state, where the order parameter r is 
computed. The procedure is repeated gradually increasing A. 

The results for r are shown in Fig.[2]against the control pa- 
rameter A for several networks characterized by different fi. 
For all values of fi, when the coupling is increased from small 
values, the incoherent solution prevails and macroscopic syn- 
chronization is not attained. This behavior persists until a cer- 
tain critical value A c (/z) is crossed. At this point some ele- 
ments lock their relative phase and synchronized nodes form. 
This constitutes the onset of synchronization. Beyond this 
value, the population of oscillators splits into a partially syn- 
chronized state contributing to r and a group of nodes whose 
natural frequencies are too spread as to be part of the coher- 
ent pack. Finally, after further increasing the value of A, more 
and more nodes get entrained around the mean phase and the 
system settles in a completely synchronized state where r « 1 
(not shown). 

A comparison between the results for different values of 
fi (and thus different (c) and (L) values) indicate several in- 
teresting features of the synchronization process. First, it is 
remarkably that when the clustering coefficient increases, the 
system reaches complete synchronization at higher values of 
the coupling. This result agrees with the results reported in 
lE7ll . where a different network model able to generate topolo- 
gies with a tunable clustering coefficient was implemented. 



At this point, one may ask whether the effects are only due 
to the influence of (c) or to the increase of the average path 
length [35] (note that the model implemented in 1E7I1 does not 
explore this possibility). Unfortunately, the two factors are 
generally linked together so they can not be considered sep- 
arately. However, as stated previously, a closer look at Fig. 
[^reveals that there is a region of the parameter \i where the 
clustering coefficient grows while the average shortest path 
length remains almost constant. This corresponds to the in- 
terval 0.4 < p, < 1.0 approximately. Going back to Fig.|2] 
the behavior of r in this interval of fx reveals that synchroniza- 
tion is almost unaffected. In fact, the r(fi) curves lie slightly 
above that corresponding to the BA limit. Therefore, though 
the above comparison is not conclusive, it seems that the de- 
layed transition to complete synchronization is mainly due to 
the effect of the increase in (L) at smaller values of fi rather 
than to the increase in (c). This conclusion is further sup- 
ported by a direct comparison of the results in Fig. |2] with 
those reported in l27ll . where the authors explored a region 
with higher values of (c) (up to 0.7) and the profile of r(A) is 
almost the same as ours. 

The second region of interest is the onset of synchroniza- 
tion. From Fig. [2] it is difficult to elucidate how the critical 
point for the BA limit compares with those at values of (j, < 1. 
At first glance, it seems that A c (/x) shifts rightward as the pa- 
rameter fx is decreased below 1. However, a more detailed 
analysis shows that it is indeed the contrary. To this end, we 
have performed a finite size scaling analysis that allows to de- 
termine the critical points A c (^). We assume a scaling relation 
of the form 

r = N- a f(NP(\-\ c )), (4) 

where f(x) is a universal scaling function bounded as x — > 
±00 and a and j3 are critical exponents to be determined. The 
estimation of A c can then be done by plotting N a r as a func- 
tion of A and tuning a for several system sizes N until the 
curves cross at a single point, the critical one. 

The results of the FSS analysis are shown in Fig. for 
different values of fi (from top to bottom and from left to 
right n = 0.05,0.15,0.50,0.60). The insets show a blow- 
up around the critical points A c (/x). Although the differences 
in the critical points at different values of fi are small, they are 
certainly distinguishable. In fact, the higher the value of fj,, the 
higher the critical point. That is, when the clustering coeffi- 
cient and the average path length grow with respect to the BA 
network, the onset of synchronization is anticipated. More- 
over, taking into account that the increase in (L) is likely to 
inhibit synchronization, one may hypothesize that the effects 
of the clustering coefficient prevail in this region of the pa- 
rameter A. To check this hypothesis, we have also included 
in Fig.|5]the analysis performed for \i = 0.50 and \i = 0.60. 
As pointed out before, for these values, the differences can 
only arise from the variations of the clustering coefficient as 
the average path length remains constant in this region of the 
parameter \i. The critical points, although very close to each 
other, are clearly different. Therefore, the main contribution 



4 



5.0 



4.0 



3.0 



s 

z 



2.0 



1.0 



0.0 



1.4 
1.1 

0.8 
0.6 
0.3 





1 /I 

f ' 




& , □□ -- 




! ^^— \=0.050~ 




• 1000 
i 2000 

• 5000 

. 10000 
*' 

A 



/ 



0.04 0.05 0.06 0.07 0.08 ♦ 

t ♦ m . • ■ ^ + 

, •»•• 



3.0 



2.0 



OS 

B 

2 



1.0 




0.04 0.05 0.06 0.07 



<x=0.220(4) 



0.04 



0.06 0.08 



1.0 



0.6 



Pi 
Z 



0.4 



0.2 



0.0 

0.10 0.02 



1.0 



WMHt!MIM ,!! "' : 



» 1000 

1 2000 i 

• 5000 a A 

l 10000 i> 

j ♦ - 



0.04 



0.06 
X 



0.08 



0.10 




0.0 — 
0.04 0.05 0.06 0.07 0/ 



a=0. 100(3) 



0.0 



< 1000 
i 2000 
■5000 
. 10000 *C*«J 

A* 

4:* 



7t 



♦ 



0.6 



o 

z 



0.02 



0.04 



0.06 
X 



0.08 



0.4 



0.2 



0.0 



0.4 
0.3 
0.2 



0.0 

0.04 0.05 0.06 0.07 0.08 



1 






' M H 




- 












— \=0.058- 


1 


1 i 1 i 



i 1000 
i 2000 
•5000 
. 10000 



A »Mf» 



a=0.090(4) 



!!! 




X* m'mr 



0.10 0.03 0.04 0.05 0.06 0.07 0.08 0.09 

X 



FIG. 3: (color online) Finite size scaling analysis for several values of /j,. From top to bottom and from left to right the values of fi are: 0.05, 
0.15, 0.50 and 0.60. In each panel, it is represented the rescaled order parameter against the control parameter A. The insets are a zoom to the 
regions around the critical points A c (/x). The data are averaged over at least 100 realizations for each value of A. The sizes of the networks, 
the critical points A c (/i) at which the onset of synchronization takes place, as well as the values of the critical exponents a are those indicated 
in the plots. See the main text for more details. 



to the onset of synchronization at low values of A comes from 
the raising of the clustering coefficient. 



DISCUSSIONS AND CONCLUSIONS 

Rounding off, our results point to a nontrivial dependence 
between the clustering coefficient and the average path length, 
and the synchronization patters of phase oscillators. Sepa- 
rately, the onset of synchronization seems to be mainly de- 
termined by (c), promoting synchronization at low values of 
the coupling strength with respect to networks not showing 
high levels of structural clustering. On the other hand, when 
the coupling is increased beyond the critical point, the effect 
of (L) dominates and the phase diagram is smoothed out (a 
sort of stretching), delaying the appearance of the fully syn- 
chronized state. These results confirm and complement those 
anticipated in Izvll and show that general statements about 
synchronizability using the MSF are misleading. Whether or 
not a system is more or less synchronizable than others show- 



ing distinct structural properties is relative to the region of the 
phase diagram in which the system operates f29 u30ll . 

In summary, we have shown that synchronizability of com- 
plex networks is dependent on the effective coupling A among 
oscillators, and on the properties of the underlying network. 
For small values of A, the incoherent solution r = first desta- 
bilizes as the clustering coefficient is higher, while the coher- 
ent solution r = 1 is promoted when both the structural clus- 
tering and the average path length are small. Finally, we point 
out that our results are also consistent if a different local order 
parameter is considered l3oll . Moreover, though these results 
have been obtained for phase oscillators, we think that they 
should hold for other nonlinear dynamical systems as well. It 
would be interesting to check this later hypothesis in future 
works. 
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